Modelos autorregresivos integrados de media móvil#
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import statsmodels.tsa.api as smtsa
import statsmodels.tsa.arima.model as arima_model
from statsmodels.tsa.arima.model import ARIMA
import warnings
warnings.filterwarnings("ignore")
url = 'https://raw.githubusercontent.com/evgomez98/wind_speed/main/wind_dataset.csv'
df = pd.read_csv(url)
from statsmodels.tsa.stattools import adfuller
from numpy import log
Prueba de Dickey-Fuller#
Iniciamos realizando la prueba de Dickey-Fuller, donde \(H_0\): La serie no es estacionaria.
result = adfuller(df['WIND'])
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
ADF Statistic: -8.413013
p-value: 0.000000
Con un \(P-Value\)= 0,00 se rechaza \(H_0\). Es decir, la serie de tiempo de la velocidad del viento es estacionaria.
Autocorrelación y Autocorrelación Parcial#
def plotds(xt, nlag=365, fig_size=(12, 10)):
if not isinstance(xt, pd.Series):
xt = pd.Series(xt)
plt.figure(figsize=fig_size)
layout = (2, 2)
ax_xt = plt.subplot2grid(layout, (0, 0), colspan=2)
ax_acf = plt.subplot2grid(layout, (1, 0))
ax_pacf = plt.subplot2grid(layout, (1, 1))
xt.plot(ax=ax_xt)
ax_xt.set_title('Time Series')
plot_acf(xt, lags=365, ax=ax_acf)
plot_pacf(xt, lags=30, ax=ax_pacf)
plt.tight_layout()
return None
plotds(df['WIND'])
El primer diagnóstico sobre la gráfica ACF se trata del patrón sinusoidal que presenta, mostrando autocorrelación y sugiriendo estacionalidad en la serie de tiempo. Por otro lado, la función de autocorrelación parcial muestra cortes significativos en múltiples rezagos lo que sugiere que varios términos de autorregresión podrían ser relevantes para el modelo. Debido a la combinación de cortes en rezagos bajos y altos, es posible que sea necesario consideras realizar pruebas de modelos ARIMA con diferentes valores de p y q, también considerar modelos SARIMA que puedan capturar tanto los patrones de autorregresión como los patrones estacionales en los datos.
División del conjunto de datos#
Para el modelamiento, se selecciona un año completo para el testeo.
tau_test = 365
tau_val = 365
train = df['WIND'][:-(tau_val + tau_test)].copy()
val = df['WIND'][-(tau_val + tau_test):-tau_test].copy()
test = df['WIND'][-tau_test:].copy()
print(f"Train: {len(train)}, Validation: {len(val)}, Test: {len(test)}")
Train: 5844, Validation: 365, Test: 365
Construcción del Modelo#
Para la selección del mejor orden de los componente del modelo, se automatiza a través de un loop que elige con base en el criterio de Akaike.
from sklearn.metrics import mean_absolute_error
best_mae = np.inf }
best_order = None
best_mdl = None
pq_rng = range(1,5)
d_rng = range(0,3)
for p in pq_rng:
for d in d_rng:
for q in pq_rng:
try:
tmp_mdl = ARIMA(train, order=(p,d,q)).fit()
val_pred = tmp_mdl.forecast(steps=len(val))
tmp_mae = mean_absolute_error(val, val_pred)
if tmp_mae < best_mae:
best_mae = tmp_mae
best_order = (p,d,q)
best_mdl = tmp_mdl
except:
continue
Cell In[8], line 3
best_mae = np.inf }
^
SyntaxError: unmatched '}'
import pickle
with open('best_arima_model.pkl', 'wb') as f:
pickle.dump((best_order, best_mdl, best_mae), f)
print('hqic: {:6.5f} | order: {}'.format(best_mae, best_order))
hqic: 4.01136 | order: (1, 0, 2)
from statsmodels.graphics.tsaplots import plot_predict
Modelo ARMA#
Entrenamos el modelo, segun el mejor orden, en este caso 1, 0, 2. Es decir, un proceso ARMA:
model = ARIMA(train, order= best_order)
model_fit = model.fit()
model_fit.plot_diagnostics(figsize=(14,10));
A través del correlograma se puede concluir independecia en los residuales.
plt.rcParams.update({'figure.figsize': (10,6)})
fig, ax = plt.subplots();
plot_predict(model_fit, 2, ax=ax);
plt.show();
Mediante la función forecast_accuracyautomatizamos el calculo de las metricas de precisión del modelo:
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy.stats import normaltest
index = 'ARMA (1, 0, 2) train'
train_result = ARIMA(train, order=best_order).fit()
train_pred = train_result.predict(start=train.index[0], end=train.index[-1])
mae = mean_absolute_error(train, train_pred)
sse = np.sum((train - train_pred) ** 2)
mape = np.mean(np.abs((train - train_pred) / train)) * 100
msd = mean_squared_error(train, train_pred)
r2 = r2_score(train, train_pred)
ljung_box = acorr_ljungbox(train - train_pred, lags=[90], return_df=True)
normality_test_stat, normality_p_value = normaltest(train - train_pred)
df_acc = pd.DataFrame({'MAE': [mae],
'SSE': [sse],
'MAPE': [mape],
'MSD': [msd],
'R2': [r2],
'Ljung-Box (p-value)': [ljung_box['lb_pvalue'].iloc[0]],
'Normalidad (p-value)': [normality_p_value]},
index= [index])
df_acc.head()
| MAE | SSE | MAPE | MSD | R2 | Ljung-Box (p-value) | Normalidad (p-value) | |
|---|---|---|---|---|---|---|---|
| ARMA (1, 0, 2) train | 3.163859 | 93215.764275 | inf | 15.950678 | 0.349119 | 2.398775e-13 | 2.992553e-66 |
index = 'ARMA (1, 0, 2) val'
val = val.reset_index(drop=True)
val_result = ARIMA(val, order=best_order).fit()
val_pred = val_result.predict(start=val.index[0], end=val.index[-1])
mae = mean_absolute_error(val, val_pred)
sse = np.sum((val-val_pred) ** 2)
mape = np.mean(np.abs((val-val_pred) / val)) * 100
msd = mean_squared_error(val, val_pred)
r2 = r2_score(val, val_pred)
ljung_box = acorr_ljungbox(val-val_pred, lags=[10], return_df=True)
normality_test_stat, normality_p_value = normaltest(val- val_pred)
df_acc = pd.DataFrame({'MAE': [mae],
'SSE': [sse],
'MAPE': [mape],
'MSD': [msd],
'R2': [r2],
'Ljung-Box (p-value)': [ljung_box['lb_pvalue'].iloc[0]],
'Normalidad (p-value)': [normality_p_value]},
index= [index])
df_acc.head()
| MAE | SSE | MAPE | MSD | R2 | Ljung-Box (p-value) | Normalidad (p-value) | |
|---|---|---|---|---|---|---|---|
| ARMA (1, 0, 2) val | 3.246131 | 5883.128472 | 36.519955 | 16.11816 | 0.304795 | 0.789485 | 0.000068 |
index = 'ARMA (1, 0, 2) test'
test = test.reset_index(drop=True)
test_result = ARIMA(test, order=best_order).fit()
test_pred = test_result.predict(start=test.index[0], end=test.index[-1])
mae = mean_absolute_error(test, test_pred)
sse = np.sum((test - test_pred) ** 2)
mape = np.mean(np.abs((test - test_pred) / test)) * 100
msd = mean_squared_error(test, test_pred)
r2 = r2_score(test,test_pred)
ljung_box = acorr_ljungbox(test - test_pred, lags=[10], return_df=True)
normality_test_stat, normality_p_value = normaltest(test - test_pred)
df_acc = pd.DataFrame({'MAE': [mae],
'SSE': [sse],
'MAPE': [mape],
'MSD': [msd],
'R2': [r2],
'Ljung-Box (p-value)': [ljung_box['lb_pvalue'].iloc[0]],
'Normalidad (p-value)': [normality_p_value]},
index= [index])
df_acc.head()
| MAE | SSE | MAPE | MSD | R2 | Ljung-Box (p-value) | Normalidad (p-value) | |
|---|---|---|---|---|---|---|---|
| ARMA (1, 0, 2) test | 3.325192 | 6441.950485 | 60.291396 | 17.649179 | 0.355096 | 0.756981 | 0.001823 |
Para el modelo ARMA(1, 0, 2) en el conjunto de prueba, el MAE de 3.33 indica que las predicciones tienen una desviación promedio de 3.33 nudos respecto a los valores observados, lo cual es aceptable en cuanto a precisión. La SSE de 6441.95 revela un error acumulado moderado, mientras que el MAPE de 60.29% muestra que la precisión es limitada en términos porcentuales. La MSD de 17.65 sugiere que el modelo produce errores significativos en sus predicciones. Un R² de 0.36 refleja que el modelo logra capturar una porción modesta de la variabilidad en la serie temporal. Además, el p-valor de la prueba de Ljung-Box (0.757) sugiere que no hay autocorrelación significativa en los residuos, favoreciendo la aleatoriedad de estos. Sin embargo, el p-valor de la prueba de normalidad (0.0018) indica que los residuos no se distribuyen normalmente, lo cual podría impactar la consistencia de los intervalos de confianza y la validez general del modelo. En conjunto, el modelo ARMA(1, 0, 2) ofrece un rendimiento aceptable, aunque con oportunidades de mejora en precisión y ajuste.
import plotly.graph_objects as go
def plot_model(train, val, test, test_pred, title):
train = np.array(train)
val = np.array(val)
test = np.array(test)
test_pred = np.array(test_pred)
train_index = np.arange(0, len(train))
val_index = np.arange(len(train), len(train) + len(val))
test_index = np.arange(len(train) + len(val), len(train) + len(val) + len(test))
test_pred_index = test_index
fig = go.Figure()
fig.add_trace(go.Scatter(x=train_index, y=train, mode='lines', name='Train', line=dict(color='blue')))
fig.add_trace(go.Scatter(x=val_index, y=val, mode='lines', name='Validation', line=dict(color='orange')))
fig.add_trace(go.Scatter(x=test_index, y=test, mode='lines', name='Test', line=dict(color='green')))
fig.add_trace(go.Scatter(x=test_pred_index, y=test_pred, mode='lines', name='Prediction', line=dict(color='red', dash='dash')))
fig.update_layout(
title=title,
xaxis_title="Observaciones",
yaxis_title="Valores",
legend=dict(x=0, y=1, traceorder="normal"),
plot_bgcolor='rgba(0,0,0,0)',
width=900, height=600
)
fig.show()
plot_model(train, val, test, test_pred, 'ARMA Model 1, 0, 2')
La función arima_rolling nos proporciona un bucle para predecir a partir de medias moviles, de modo que el modelo se adapte a la evolución de los datos:
def arima_rolling(train, val, test, best_order):
history = list(train) + list(val)
predictions = []
for t in range(len(test)):
model = ARIMA(history, order=best_order)
model_fit = model.fit()
output = model_fit.forecast()
yhat = output[0]
predictions.append(yhat)
obs = test[t]
history.append(obs)
print(f'Predicted={yhat:.3f}, Expected={obs:.3f}')
return predictions
Ahora pronosticamos y comparamos con los datos de prueba:
test_pred = arima_rolling(train, val, test, best_order)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[1], line 1
----> 1 test_pred = arima_rolling(train, val, test, best_order)
NameError: name 'arima_rolling' is not defined
index = ARMA rolling (1, 0, 2) test
mae = mean_absolute_error(test, test_pred)
sse = np.sum((test - test_pred) ** 2)
mape = np.mean(np.abs((test - test_pred) / test)) * 100
msd = mean_squared_error(test, test_pred)
r2 = r2_score(test,test_pred)
ljung_box = acorr_ljungbox(test - test_pred, lags=[10], return_df=True)
normality_test_stat, normality_p_value = normaltest(test - test_pred)
df_acc = pd.DataFrame({'MAE': [mae],
'SSE': [sse],
'MAPE': [mape],
'MSD': [msd],
'R2': [r2],
'Ljung-Box (p-value)': [ljung_box['lb_pvalue'].iloc[0]],
'Normalidad (p-value)': [normality_p_value]},
index= [index])
df_acc.head()
| MAE | SSE | MAPE | MSD | R2 | Ljung-Box (p-value) | Normalidad (p-value) | |
|---|---|---|---|---|---|---|---|
| ARMA (1, 0, 2) test | 3.326124 | 6425.376786 | 61.322275 | 17.603772 | 0.356756 | 0.725282 | 0.001373 |
Para el modelo ARMA(1, 0, 2) en el conjunto de prueba, el valor de MAE de 3.33 sugiere una desviación promedio de aproximadamente 3.33 nudos en las predicciones, mostrando una precisión aceptable en términos absolutos. La SSE, con un valor de 6425.38, denota un error acumulado moderado en el ajuste, mientras que el MAPE de 61.32% indica que el modelo no es muy preciso en términos relativos o porcentuales. La MSD de 17.60 confirma una magnitud significativa en los errores, y el coeficiente de determinación \(R^2\) de 0.36 muestra que el modelo logra capturar parte de la variabilidad en los datos, aunque con margen de mejora. El p-valor de la prueba de Ljung-Box (0.725) sugiere una ausencia de autocorrelación significativa en los residuos, favoreciendo la aleatoriedad en estos. Sin embargo, el p-valor de la prueba de normalidad (0.0014) sugiere que los residuos no siguen una distribución normal, lo cual podría afectar la robustez de los intervalos de confianza. En conjunto, el modelo ofrece un rendimiento razonable, aunque podría mejorarse en términos de precisión relativa y ajuste de los residuos.
plot_model(train, val, test, test_pred, 'ARMA rolling Model 1, 0, 2')
pred_error = np.array(test) - np.array(test_pred)
from pandas.plotting import autocorrelation_plot
fig = plt.figure(figsize=(5.5, 5.5))
autocorrelation_plot(pred_error, color='b')
plt.title('ACF for Prediction Error');
Tal como el Ljung-Box no rechaza la independenicia de los residuales del modelo ARMA(1,2), se puede apreciar en la gráfica este comportamiento.
Modelo SARIMA#
import statsmodels.api as sm
best_aic = float('inf')
best_pdq = None
best_PDQ = None
best_mdl = None
for p in range(5):
for d in range(3):
for q in range(5):
try:
tmp_mdl = sm.tsa.SARIMAX(train_data, order=(p,d,q), seasonal_order=(p,d,q,365)).fit()
tmp_aic = tmp_mdl.aic
if tmp_aic < best_aic:
best_aic = tmp_aic
best_pdq = (p, d, q)
best_PDQ = (p, d, q)
best_mdl = tmp_mdl
except Exception as e:
continue
#if best_pdq and best_PDQ and best_mdl:
#with open('best_sarima_model.pkl', 'wb') as f:
#pickle.dump((best_pdq, best_PDQ, best_mdl, best_aic), f)
#print(f'aic: {best_aic:6.5f} | order: {best_pdq} | seasonal_order: {best_PDQ}')
#else:
#print("No se encontró ningún modelo adecuado.")
El loop para el modelo SARIMA no se ejecuta por el riesgo de sobrecalentamiento del equipo.